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Abstract: We consider the modeling of the dynamics of the chemostat at its very source. The 
chemostat is classically represented as a system of ordinary differential equations. Our goal is to 
establish a stochastic model that is valid at the scale immediately preceding the one corresponding 
to the deterministic model. At a microscopic scale we present a pure jump stochastic model that 
gives rise, at the macroscopic scale, to the ordinary differential equation model. At an intermediate 
scale, an approximation diffusion allows us to propose a model in the form of a system of stochastic 
differential equations. We expound the mechanism to switch from one model to another, together 
with the associated simulation procedures. We also describe the domain of validity of the different 
models. 
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Modeles stochastiques du chemostat 



Resume : Nous reprenons la modelisation de la dynamique du chemostat a sa source. Le chemo- 
stat est classiquement represente par un systeme d'equations differentielles. Notre objectif est 
d'etablir un modele stochastique qui est valable a l'echelle qui precede immediatement celle qui 
correspond au modele deterministe. Partant d'une echelle microscopique, nous presentons un 
modele stochastique de sauts purs qui conduit, a l'echelle macroscopique, au modele d'equation 
differentielle. A une echelle intermediaire, une approximation diffusion nous permet de proposer 
un modele sous la forme d'un systeme d'equations differentielles stochastiques. Nous detaillons les 
techniques qui permettent de passer d'une echelle a une autre ainsi que de simuler ces differents 
modeles. Nous decrivons egalement les domaines de validite des differents modeles. 

Mots-cles : equations differentielles stochastiques, chemostat, processus de saut, approximation 
diffusion, methode "tau-leap", methode de Monte Carlo, algorithme de Gillespie 
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1 Introduction 

The evolution of the state of a single species/single substrate chemostat is usually described by a 
set of ordinary differential equations (ODE) derived from a mass balance principle, see [25]. More 
precisely, if s(t) denotes the concentration of nutrient (substrate) and b(t) the concentration of the 
organism (biomass) at time i(expressed in g/L), then the couple x(t) — (b(t), s(t))* is the solution 
of the following ODE [25]: 

b(t) = [v(s(t))-D}b(t), (la) 
s(t) = -kfi(a(t)) b(t) + D [s ia - s(t)} (lb) 

where D > is the dilution rate, s ,n > the substrate concentration in the influent, and k > the 
stoichiometric coefficient. The initial condition lies in the positive orthant, that is 6(0) > and 
s(0) > 0. Equation (1) will also be denoted: 

x(t) = f(x(t)) . 

The specific growth rate function /z(s) is non-negative; we suppose that /i(0) = 0, /i(s) > for 
s > 0, fi(s) < /i max < oo and that it is continuous at 0. Commonly used models are the Monod 
model (uninhibited growth) and the Haldane model (inhibited growth) that reads respectively: 

s s 

ll(s) = u max , u(s) = u max T ■ (2) 

W K + s PW Pmax fc s + s + % 

This approach relies on the fact that the stochastic effects can be neglected, thanks to the law of 
large numbers, or at least can be averaged out. Although this level of description is sufficient for a 
number of applications of interest, it could be a valuable way of accounting for the stochastic nature 
of the system. Indeed, at small population sizes the chemostat could present stochastic behaviors, 
also the accumulation of small perturbations in the context of multi-species could not be neglected. 
Moreover, whereas the experimental results observed in well mastered laboratory conditions match 
closely the ODE theoretical behavior, a noticeable difference may occur in operational conditions. 
In these cases, stochastic features may not be neglected. We aim to build a model that still relies 
on a mass balance principle and that encompasses the useful stochastic information. 

Many works [26, 7, 15] propose to superpose a stochastic term on Equation (1) in order to 
model the uncertainty on the phenomenon, principally due to imprecise experimental conditions. 
Paradoxically, this amounts to the addition of an ad hoc perturbation to a model that has been 
obtained by neglecting these perturbations. We propose instead to consider the stochastic aspect 
at the very beginning of the modeling process, and to determine the conditions under which it 
is insignificant. This approach is not individual-based per se, as it starts from the macroscopic 
model (1). However, the first stochastic model proposed will be described at the individual level. 
This method will allow for a justification of the specific structure of the stochastic perturbation 
that affects the mean behavior. More generally, we will outline a modeling strategy based on many 
available tools, either stochastic or deterministic, depending on the regularity of the phenomenon 
to be modeled. In this paper we focus on the modeling and simulation process rather than on the 
mathematical developments; moreover we make use of known mathematical results. Our goal is to 
establish a stochastic model that is valid at the scale immediately preceding the one corresponding 
to the deterministic model (1). 

The paper is organized as follows: in Section 2, we recall the origin of model (1) and the 
assumptions ensuring its validity. We show that since different timescales naturally appear in the 
problem, these assumptions need to be checked at each scale. Section 3 is devoted to the different 
models: the pure jump description that will be considered as the reference model is introduced 
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in Section 3.1; the discrete time approximation, Poisson and normal, are presented in Section 3.2; 
the discrete-time normal approximation appeared to be a time discretization of a diffusion process 
given by a stochastic differential equation presented in Section 3.3. In Section 3.4 we describe 
the asymptotic results that bridge these different models. Section 4 is devoted to the associated 
simulation algorithm, Section 5 to numerical tests. 

2 Scale and geometry issues in ODE model 

An individual-based model should keep track of the position in space of each cell, together with 
their current biological states, it should also account for discrete events such as the division of a 
cell. Such a description of the system at the finest level could be of interest but unnecessary in 
view of our goals, namely to set a macroscopic model that gives account for stochastic phenomena. 
At this scale, the system is reduced to a R 2 -vector and its dynamics. 

Model (1) is obtained according to the classical approach, by choosing a small time interval At 
on which a mass balance principle is applied to the state. However, At should be large enough as 
we do not describe the dynamic at the timescale of jumps of one unit of substrate or bacteria but 
rather at the timescale of jumps of packet units. Such an interval could be called macroscopically 
infinitesimal [8]. 

Mass balance 

Let (Bt,St) denote the true concentrations at time t, assumed to be constant throughout the 
medium. The mass balance on interval [t, t + At) reads 

Bt+At -Bt = AB?° + ABr , (3a) 
St+At -St = AS ] r + AS' t " + AS° M (3b) 

where 

• A2? bio and A£>j ut are the increments of biomass due to natural growth and to the outflow 
respectively, within [t, t + At), 

• A<S f bio , A6>j" and AS™* are the increments of substrate due to the consumption by the 
biomass, the inflow and the outflow respectively, within [t, t + At). 

Since we want to obtain an ODE, we now assume that the stochastic fluctuations are negligible 
relative to the increments. Again this requires At to be large enough, so that sufficiently many 
discrete events have occurred. Moreover, At should be taken even larger in case of inhomogeneity 
of the dynamics. 

We denote by (6(t),s(t)) for t = 0, At, 2 At,... the deterministic sequence constructed by 
using the mean increments of (B t ,S t ): 

b(t + At) - b(t) = E[A^'° + A£° ut ] , 

s(t + At) - s(t) = E[A«S t bl ° + AS l t n + AS° ut ] ■ 

Next, using the mass action law for the biomass we have 

E[AB b '°] ~ n(s(t)) b(t) At , E[A5 4 bio ] ~ -k n(S(t)) b(t) At (4) 

where /x(s) is the specific growth rate and k > the stoichiometric coefficient. Note that we again 
require At to be large enough, since fj,(s) and k make sense only for a sufficiently large population 
of bacteria. Now, since we have assumed perfect homogeneity of the medium, we get: 

E[AB° t ut ] ^ -D b(t) At , E[AS; n ] ~ —D s'" At , E[A5 f out ] ^ -D s(t) At . (5) 
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Note that (4) and (5) are approximations because we have used a constant value for b(t) and s(t) 
within [t,t + At). For this approximations to be correct, none of the quantities involved should 
vary significantly within [t, t + At). We finally obtain the construction of the sequence (b(t), s(t)) 

by 

b(t + At) - b{t) = [fi(s(t)) - D] b(t) At, (6a) 
s(t + At) - s{t) = {-k fi(s(t)) b(t) + D [s in - s(t)]) At, (6b) 

Model (1) is obtained by letting At — > in System (6). However, since At is bounded from below, 
some care should be taken when this limit is achieved. System (6) can be understood as the 
discretization of (1) using an explicit Euler scheme with time-step At. Whenever there exists At 
sufficiently small, the deterministic sequence (b(t), s(t)) will be close to model (1), sampled at time 
0,At, 2 At,... 

Geometry and scales 

The mass balance established in (3) features five terms that can be gathered according to the 
three sources of variations. This gives rise to a geometric structure that can be emphasized by 
writing (1) under the form: 

£ ( S ) - ti'WW) i±u) +D (&) -D C$) (7) 

biology inflow outflow 

However, whereas the geometry is well captured, the timescale of the five original terms is not 
readable in (1) nor (7). Indeed, the fact that the approximations in (4) and (5) may be of different 
quality for each individual term is not exploited at all. 



3 Models at different scales 

In the previous section, we mentioned that the lower bound for At is related to the size of the 
population and to the regularity of the phenomenon. Often, the experimental conditions are such 
that this bound is low enough, so that System (6) is correctly approximated by (1) sampled with 
period At. If a smaller period is to be considered, then the conditions under which (6) has been 
obtained are not fulfilled. Particularly, the stochastic fluctuations should be accounted for. 

We now introduce a stochastic process built on the same premise, that is a mean mass balance 
principle at a given At. This model will have (f) as a fluid limit as At goes to 0. This latter 
model suitably features the geometry of the chemostat but, as a limit model, cannot feature all 
its natural scales. The proposed stochastic models will respect both the geometry and the natural 
scales of the chemostat. We first establish a pure jump process representation of the chemostat 
at a microscopic scale, then we derive a diffusion process representation which will be valid at 
mesoscopic and macroscopic scales. 

3.1 Pure jump model X t = (B t , S t )* 

Even if do not aim at deriving an individual-based model, we try to preserve the discrete feature 
in the dynamics. We achieve this by considering only aggregated jumps obtained by adding up 
small and frequent jumps resulting from individual events. The resulting stochastic process will 
be a pure jump process X t — (B t , St)* , fully determined by its jumps and the corresponding jump 
rates; the state variable will be denoted x — (b, s)* . 
In view of (3), we are led to consider five jumps: 
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® biology term: biomass increase of size v\{x) at rate Ai(a:); 

© biology term: substrate decrease of size 1*2(2;) at rate A2(x); 

© inflow term: substrate inflow of size 1*3 (x) at rate A3 (a;); 

@ outflow term: biomass outflow of size 2/4(2;) at rate X^x); 

© outflow term: substrate outflow of size v§ (x) at rate \^,(x); 

(see Figure f ). It remains to set the jump size rates so as to comply with the mass balance principle 
and the stochastic mass action law. 

For a macroscopically infinitesimal At, denote by AX t b bl °, AX t ! ' bio , AX t 3,in , AX t b ' out , AX t a ' bio the 
cumulated jump of type ©, ©, ©, ©, © respectively, on state process X t within the time interval 
[t,t + At). 

We first focus on the first two expressions. The stochastic mass action law [28] requires 

E[AXl bi °\X t =x]~ (mM* At) 5 

E[AXt ia \X t = x}~(_ kK % bAt ) . 

Now notice that, for small At, the number of jumps of type © (resp. ©) within [t,t + At) is 
approximately V(X\ (x) At) (resp. V(X 2 (x) At)), so that 

E[AXp bIO \X t = x] ~ Ai(x) Atvi(x), 
E[AX s t bio \X t = x] ~ \ 2 (x) Atv 2 (x) . 

So we are looking for (A^(x), Vi(x)) satisfying: 

X 1 (x)u 1 (x) = (^ b ) and A a (x)^(x) = (_ fcp ° W6 ) . (8) 

We therefore introduce the scale parameters K\ and K 2 and we choose: 

Ai(x) = ifi/x(s)6, Ul =(^T^, 

A 2 (x) = K 2 k fj,(s) b u 2 = -(^-J. 

this choice is not unique and will be explain later in Section 3.5. 

Here by "scale" we mean that jumps due to © will be of magnitude and the corresponding 
rates will be of magnitude Ki. Large Ki yields frequent and small jumps. Using the Poisson argu- 
ment mentioned above, we see that these scale parameters Ki do not act on the mean values of the 
increments but on their variances (large Ki will correspond to small variances). The Ki's can thus 
be regarded as tuning parameters quantifying the uncertainty or regularity of the corresponding 
source of variation. 

Reproducing this discussion with the three other types of jumps, and considering only admis- 
sible jumps (in the positive orthant), we obtain a pure jump Markov process with rate coefficients 
Xi(x) and associated jumps Vi(x) defined in Table 1. 

About scales parameters 

Let m b and m B denote the representative masses of a single bacteria and of a single molecule of 
substrate. Typically m b 3> m s (e.g. m b ~ I0 6 m 3 ). Hence: 

< Ki <—, for i = 1,4 and < K % < — , for i = 2, 3, 5 . 
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Table 1: Rates and jumps of the five basic mechanisms of the pure jump process. Note that 
the jumps Vi(x) essentially do not depend on x except for the negative jumps near the border 
{x = (6, s)* G R^; b = or s = 0}. 
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Figure 1: In this model, from a position x = (b,s)* the process could jump according to 5 mech- 
anisms (2 due to the biology, 1 inflow, and 2 outflows), the basic jump © has a length for 
i = 1,...,5. 



In most cases Ki <C Kj for i = 1, 4 and j = 2, 3, 5, but it is possible to adjust the coefficients 
i^'s to the specific application considered. For example K2 will be large in laboratory experi- 
mental conditions, but for a real implementation the substrate inflow concentration could have a 
large variance. Also for the outflow, in regular conditions K4 and K$ could be large, but in bad 
mixing conditions they could be smaller. Finally K\ could be smaller than K4, as the biomass 
concentration increase presents more variance than the substrate decrease (which is more regular 
as it is related to the diffusion of substrate across cell membranes). 

As proved later in Lemma A.l, the jumps Vi(x) are essentially constant and equal to: 

-=ff), *«-(*)• **-(*)■ •**-(£)• <« 



Representation of X t 

The constructive description of the process X t that has been just presented would be used for 
simulation purposes, see Section 4.1. Nevertheless, it should be completed by a more comprehen- 
sible and synthetic representation. This will require some mathematical developments which we 
summarize now and that are detailed in Appendix A. 
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First we should notice that the jump process X t can be represented as the following (jump) 
SDE: 

X t = X + Y / f Vi{X u -)l {v < Xi{x )} N\duxdv) (10) 

i=1 «/(0,t]x[0,oo) 

where N l are independent random Poisson measures with intensity measure du x dv (Lebesgue 
measure) . 

The process X t can be described as a Markov process with infinitesimal generator A<f>(x) — 
lim^o 7 — (f>(x)] where Xf is the process X t starting from x. This operator completely 

characterizes the law of the process X t . In Appendix A we prove that this process is non explosive, 
i.e. it is defined for all t > 0; that it admits moments as soon as Xo does; and that it is solution 
of (10) see Proposition A. 3. 

Still representation (10) is opaque. We can establish that the process X t essentially admits the 
following representation: 

dB t = ^(S t )B t -DB t )dt+^ + ^0 : , (11a) 

dS t = ( -k^St) B t + D (,s'» - St)) dt + ^§ + *g + *g (lib) 

where fh\ are independent square integrable martingales with zero mean. The exact representation 
(41) differs from (11) only through terms (1 A s) and (1 A Ki b) that are equal to 1 except on a 
very limited neighborhood of the axes. 

The martingales ml are of mean and they are explicitly known, see (40), as well as their 
quadratic variation, see (42). From Equation (11), the deterministic part of the process X tl its 
drift coefficient, appears to be essentially the classical ODE (1); and the stochastic part of this 
dynamics, the martingale terms, are of order l/^/lQ. 

3.2 Discrete time approximations 
Poisson approximation X tn = (Bt n ,St„)* 

For any small At > given, let t n = n At. We propose a discrete time Poisson approximation 
(X tri ) n >o of (A t ) t >o: on the interval [t n ,t n+ i) we froze the rate functions Xi(X t ) to Xi{X tn ) so 
that we get a Poisson distribution. The jumps Vi{X t ) are also frozen to Vi{X tn ). Let X — X , 
the approximation is defined by: 

5 

x tn+1 = x tn +Y,<x tn )K(^HK)) (12) 

where ('P I l l ( j o)) ne N.i=i---5 ar e independent Poisson variables with intensities p. 
We have: 

5 

E[X tn+1 \X tn =x]=x + ^v i (x)W(P i n {At\ i {X tn ))\X tn =x] 

i=l 

5 

= x + At^ j v l (x)\ l (x) (13) 

i=l 

and let 

5 

MaO^X^aOM*). (14) 

i=l 
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Jk(x) is "essentially" the r.h.s. function f(x) of the O.D.E. (1), more precisely fk{%) — f(x) except 
near the axes (see. Lemma A.l). In other words, the infinitesimal increments of the conditional 
mean follow the O.D.E. (1). 
Also: 

cov[X tn+1 \X tn =x]=J2 covh(x) V\(At Xi(X tn ))\X tn = x} = 9\ (i 5a ) 

i=i ^ 2 ' 

with 

t\ = ^cov[P^{AtX 1 (x))} + ^(lAK 4 b) 2 cov[V^(AtX 4 (x))] 

= At{^ f i( S )b+^-JlAK i b) 2 Db}, (15b) 

t 2 = -^(lAK 2S ) 2 cov[V 2 (AtX 2 (x))] + -^cov[VUAtX 3 (x))} 

+ ^(lAK 5 s) 2 cov[V^AtX 5 (x))} 

- At { ±- 2 (1 A K 2 sf k M (s) b + ± D + ^ (1 A K 5 s) 2 D s} . (15c) 
Diffusion approximation £ tn = (/3t„,<7t n )* 

In (12), the variable V^AtX^x)) is Poisson distributed with parameter AtXi(x). When this 
parameter is large (greater than 10 or 20) then this last distribution is very close to the normal 
distribution of mean AtXi{x) and variance AtXi{x). Hence, we get a (discrete time) normal 
approximation (£t„)n>o of (X t )t>o by letting £o = Xo and, conditionally on ^t n _ 1 = x: 

5 
i=l 

where Af„ are 5 independent Gaussian random variables : 

N l n ~N(Xi{x) At, Xi(x) At) 

So conditionally on £t n = x, £t„ +1 i s normal with mean (13) and covariance matrix (15). 
Let | tn = (f3 tn ,a tn )*, given (3 tn = b and a tn = s: 

Pt n+1 = b + [/,(«) - (1 A K, b) D]bAt+ ^At^ wi + /Ai |M f flt < (16a) 

a tn+1 = s + [ - (1 A K 2 s) k /Lt(s) b + D s ia - (1 A K 5 s) D s] At 

+ y/At wl + y/At 2£wl + y/At QL^gH w l (16b) 

where w l n are i.i.d. A/"(0, 1) random variables. 

Boundary conditions 

In both approximations (12) and (16), no mechanism prevents the processes X tn or £ tn from staying 
within the positive orthant R+. An ad hoc solution is to set the concentration to whenever it 
becomes negative, see Section 4. 
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3.3 Diffusion model & = (f3 t ,a t )* 
A stochastic differential equation model 

System (16) is the Euler-Maruyama time discretization of the diffusion process £ t = (/3 t ,<r t )* 
solution of the following SDE: 

dA = [K<rt) ~ (1 A K 4 (3 t ) D] ft dt + dW i + ^KMEEE dW f 

da t = [- (I A K 2 o- t ) k fj,(a t ) /3 t + D s in - (1 A K 5 a t ) D a t ] dt 

+ ^kI^EEEAK dW 2 + JD£_ dW 3 + ^ (lAK 5 a t )i D *; ^5 

where W t 8 are independent standard Wiener processes. Note that this result can be obtained 
directly from the process X t without the help of the discrete-time approximation. Indeed the 
infinitesimal generator of process X t given by (26) is a difference operator, and by Taylor develop- 
ment, it can be approximated by a second order differential operator corresponding to a diffusion 
process [6]. 

For small K^s the last system is equivalent considering: 

dA = [»(* t ) p t -D ft] dt + y^F dw ' + \/W AW ' 

da t =[-k n{<7 t ) ft + -D s in - D a t ] dt 

then we can group the Brownian motions in the following way: 

dft = [^(a t ) - D] ft dt + yJlitsgei + °& dW? (17a) 

da t =[-k M ( CTt ) ft + £> ( S '» - a t )] dt + + n wT + ! T^ dW i ( 17b ) 

where W h and W s are independent standard Wiener processes. 

Behavior of the system of SDE's near the axes 

System (16) is the Euler-Maruyama time discretization of the SDE (17) (for large Ki's). Even if 
the diffusion approximation is only valid for large values of the biomass and the substrate, we can 
study the behavior of (17) near the axes. 

As for the discrete-time normal approximation, we should clarify the boundary conditions. As 
we well see, the component ft given by (17a) will remain positive, but the component er t given by 
(17b) could become negative. We must first require that fi(s) = for s < 0. Then, note that each 
equation of (17) is related to the following well-known CIR model for interest rates: 

Remark 3.1 (Cox-Ingersoll-Ross model) Consider the one-dimensional SDE: 

dX t = {a + bX t )dt + ayfx' t dW t , X o = x Q >0. (18) 

with a > 0, b € R, a > 0. According to [19, Prop. 6.2.4], f or a ^ Xq > 0, X is a continuous process 
taking values in R + ; and let t = inf{i > 0, X t = 0} ; then: 

(i) If a > cr 2 /2, then r — oo P x -a.s.; 

(ii) if < a < <7 2 /2 and b < then r < oo V x ~a.s.; 
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(in) ifO<a< a 2 /2 and b > then P x (t < oo) £ (0, 1). 

In the first case, X never reaches 0. In the second case X a.s. reaches the state 0, in the third 
case it may reach 0. If a = then the state is absorbing. 

It is clear that j3 = is an absorbing state for (17a), and when (3 = 0, (17b) reduces to 

dat = D ( S '» -a t )dt+ v^pT^f dWi 

and from Remark 3.1, the solution of this SDE will stay on the half-line [— s"\oo) and: 
(i) if s ia (-^ + -jj^) > jj^F then at never reaches — J|| s ,n ; 

(zi) if s ln (-^ + -^) < ^-^r then cr t reaches — s m in finite time and is reflected. 

Indeed, it is enough to apply Ito formula to at = — h and to use the Remark 3.1. Note 
that, as K$ is large, condition (z) is more realistic than condition (ii). 

To extend the definition of (17) for negative value of a, let suppose that fi(a) — for a < 0. 
As we seen, j3 t will stay non- negative and /3 = is an absorbing state. Also a t > —j^ s ia and 
for large K§ this state will be repulsive. Note that for small values of cr t , as the Ki are large, the 
diffusion term in (17b) will be small and the drift part will be dominated by Ds in so that at will 
increase fast and its probability to be negative will be small. 

The fact that the substrate concentration could be "negative" is due to the normal approxima- 
tion. This approximation is valid for large values of concentration and the validity of the diffusion 
system (17) is questionable for small concentration. Nonetheless we can study its properties. 

A possibility to get an SDE with positive solution is to consider an SDE with boundary condition 
[14, § IV- 7] by adding a local time in {a — 0} to the Equation (17b). This solution is rather artificial 
and will not be retained. 

So the solution of the system (17) remains in the domain T> = [0, oo) x [— s in , oo). The proof 
that this system admits a strong solution with pathwise uniqueness is presented in Appendix B. 

3.4 Asymptotic analysis 

The convergence of the pure jump model (10) or of the diffusion approximation (17) to the deter- 
ministic model (1) as all the Ki — > oo can be rigorously established. 

Let X^ be the pure jump model defined at the beginning of Section 3.1, or as the solution 
of the Equation (10) for a given K = f (K\, K 2 , K 3 , K^, K 5 ). Let be the solution of the SDE 
(17). Let x(t) be the EDO model solution of Equation (1). Then X{~ converges toward x(t) in the 
following way: for all T > and all S > 0, 

p( sup \\Xf -x(t)\\ >S) — >0 (19) 

as Ki — > oo for all i = 1 ■ ■ ■ 5. This result is not surprising if we consider the representation (41) of 
(Xt)t>o; it was obtained in a context of martingale convergence theorems in [17, 18] or in a more 
general context of convergence of sequences of infinitesimal generators in [6] . 

We can also prove the same type of convergence for the process £^ . Indeed, in Equation (17) 
the scale coefficients appears as \j\fK, in the diffusion part of the SDE, and the convergence 
clearly holds as all the Ki tends to infinity. 
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Figure 2: In this simplified model, from a position x — (6, s)* i/ie process could jump according to 
3 mechanisms (biology, inflow, and outflow), the basic jump ©' has a length for i = 1, . . . , 3 

3.5 Other models 

As already noticed, the choice of (Xi(x), Vi(x)) satisfying (8) is not unique. We choose not to make 
the jump sizes depend on the state value x (except for the boundary conditions), only the jump 
rates depend on x. Another possibility is to choose jump sizes that depend on the state value x. 
For example instead of the choice of Table 1, we can choose: 

Ai(x) = ifi//(s), A 2 (as) =K 2 k^{s), \ 3 (x) = K 3 D, \ 4 (x) = K 4 D , \ 5 (x) = K 5 D 

and 




(if we neglect the boundary condition). Then in place of (17) we have the following set of equations: 

dA = [ufa) fit- D [3 t ] dt + sf^-flt dWt 1 + sffjt dW? (20a) 
da t = [ - k tx{a t ) fit + D s ln - D a t ) dt 

+ -^W^ dW t + V^ s '" dW ? + \[% a * dw " ( 20b ) 

where W l are independent standard Wiener processes. 
A three components model 

Instead of the five components © to ©, we can consider a case with three independent sources of 
jump variation. This example strictly preserves the geometry (7) by considering three independent 
sources of jump variation: 

©' biology term: biomass increase and substrate decrease at scale K[; 

©' inflow term: substrate inflow at scale K' 2 \ 

©' outflow term: biomass and substrate outflow at scale 



RR n" 7458 



Stochastic models of the chemostat 



15 





©' ©' ©' 


biology inflow outflow 


rate X'i(x) 


K[ a(s) b K' 2 D K'^D 


jump v[{x) 


( $ \ / o \ /-^A 



Table 2: Rates and jumps of an ad hoc choice for three mechanisms of the pure jump process. Note 
that the third jump v' 3 (x) is radial. 

see Figure 2. 

Again the jump sizes and rates should be chosen so as to satisfy the mass balance principle and 
the stochastic mass action, with no canonical choice. An ad hoc choice is given in the following 
table: 

These jumps are now essentially equal to 

, a* 1 / 1 \ , dot- 1 / \ , dot 1 f-b\ 

This setting forces the jumps to be directed along the corresponding vector field, which is a strong 
constraint. In particular, the stoichiometry is strictly respected: the production of 1 unit of 
biomass requires exactly k units of substrate. Moreover the outflow jump is always radial, so that 
the increments of biomass and substrate are again strongly linked. Notice that for this particular 
choice of A3 and f 3 , the jump rate is constant but the jump size is not. In other words, the jump 
carries information both in the direction and the intensity of the variation. This will affect the 
qualitative behavior of the process and of its diffusion approximation, regarding extinction for 
example. 

As for our canonical model, we obtain a SDE for the diffusion approximation of the jump 
process : 

dft = [fi(a t ) - (1 A K' 3 ft) D] f3 t dt + dWi + Jgr(l A X 3 II6II) ft dW? 

da t = [ - (fc A K[ a t ) fi(a t ) Pt + D s' n - (1 A K' 3 ||&||) Dtr t ] dt 

+ V^IF^ A K 'x dW t + ^ s ' n + yf%(l A K> liej) a t dW t 5 . 

Notice that since W 1 and W 3 affect both components of £ t , the quadratic variation process (£)t 
will not be a diagonal matrix. In comparison with (7), we write the vector form of the SDE for 
large Ki's: 




biology inflow outflow 




biology inflow outflow 

The diffusion term appears as the conjunction of three perturbations acting along the three vector 
fields determined by the sources of variation. Moreover, the intensity of the noise could be different 
for each type of perturbation. Considering this model could therefore be of interest, if the geometric 
interpretation of the noise is meaningful, see [16]. 
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Comparison with the Imhof-Walcher model [15] 

We finally mention that the diffusion model appearing in [15], is obtained from (20) by letting 
K\ = if 2 = K$ = which leads to: 



The choice of these coefficients is justified in [15] by constructing an approximating Markov chain, 
and then taking the limit as the sampling rate goes to 0. This model will be compared to the 
diffusion approximation (17) model on a simulation test in Section 5.4. 

4 Simulation algorithms 

We presented several models for the chemostat system: the pure jump model (X t )t>o could be 
considered as a detailed model at the microscopic scale. The Poisson approximation (Xt n ) n eN given 
by (12) and the normal approximation (£t„)«.£N given (16) are constant time step approximation of 
the pure jump process. Finally the diffusion process (£t)t>o solution of the SDE (17) is a continuous 
time approximation of the pure jump process. 

The now present the three associated simulation algorithms that will be valid at different scales. 

4.1 Pure jump model 

The pure jump model in continuous time described in Section 3.1 can be exactly simulated thanks 
to the Gillespie algorithm, also called stochastic simulation algorithm, described in Algorithm 1. 

When the rate coefficients Xi(x) are large the time increment will be small and the Gillespie 
algorithm is impractical. As the scale coefficients Ki are large, the Xi(x), i 7^ 3, are large only 
when j3 and a are small; A3 (a;) will remain large as it does not depend on x. 

4.2 Poisson approximation 

The simulation of the previous model could be cumbersome for very high rates of event. In this 
case it is desirable to use the fixed time step Poisson approximation method (12) also called tau- 
leap [9] . Recently many papers have addressed the numerical analysis of this approximation scheme 
[23, 20, 1]. In this method the time step should be small enough so that it fulfills the following 
"leap condition": the state change in any leap should be small enough that no rate function Xi{x) 
will experience a macroscopically significant change in its value, that is: 



for i = 1 • • • 5, where < e <C 1 is an error control parameter. 

For this method to be practicable [10] proposed an automatic and simple way of determining 
the largest time step At compatible with the leap condition. Define: 



Ap t = - D] ft df + & P t dW? 

do-t = [ - k ii(a t ) Pt + D (s' n - a t )] dt + c s a t dW t B 



(21b) 



(21a) 



\i(x + J2i> "i'WPni^t \i>(x))) - \i(x) <eX{x) 



(22) 



5 




(23a) 



i'=i 



5 




(23b) 



i'=i 
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t <— 0, x <— xq 

save (t, x) 

while t < T max do 

compute Xi(x) °/« see Table 1 

A = Ei=i ^i(ar) 

At ~ Exp(A) 7, exponential distribution 

u ~ U[0, 1] °/o uniform distribution 

* <- t + At 

if u < Ai(x)/A then 

x <— x + v\ % biomass reproduction 
else if u < {Xx(x) + \2{x)}/A then 

x ■h- [x — + °/» consumption 
else if u < {Xi(x) + X 2 (x) + A3(x)}/A then 

x <— x + P3 7« substrate inflow 
else if it < {Xi(x) + X 2 (x) + Aa(^) + A4(x)}/A then 

a; <— [x — ^4]+ 7. biomass outflow 
else 

x <— [x — ^5]+ 7. substrate outflow 
end if 
save (t, x) 
end while 

Here ^ = ( V* ) , p 2 = ( ^ ) , ^ = ( ^ ) , *7 4 = ( V*« ) , P 5 = ( ^ ) and [x] + is 
the projection on the positive quadrant: [x] + = [(£)], = (^vo)' 

Algorithm 1: Gillespie algorithm (or stochastic simulation algorithm). 



t <— 0, a; <— xo 

save (t, x) 

while f < T max do 

compute Xi(x) 7, see Table 1 

A = Ei=i A iO) 

compute m,i(x), Vi(x) °L see (23) 

At <- min i= i... 5 {e X/\rrn(x)\ , e 2 A 2 /vi(x)} 

t <- 1 + At 

~ Poisson(Ai(x) At) for i = 1 • • • 5 

2; «- [x + Pi Pi - ^2 + P3 V-i - Vi Va - v 5 7> 5 ] + 

save (t, x) 
end while 

Here v x = ( V* ) , „ 2 = ( ^ ) , p 3 = ( ^ ) , p 4 = ( V*« ), ^ = ( 1/ ^ ) and [x} + is 
the projection on the positive quadrant: [x]+ = [(£)]. = (fvo)' 

Algorithm 2: Poisson approximation or tau-leap method. 
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for i, i' = 1 • • • 5, and let 

a1 =,s'° 5 {kUi' Mir; (24 » 

where e is an error control parameter (0 < £ <C 1), see Algorithm 2. Note that in the original 
context the jumps Vi(x) do not depends on x, but in our situation they do not essentially depend 
on x, the dependence on x was introduced to handle the jump near the axes in order to avoid 
negative concentration. 

4.3 Diffusion (normal) approximation 



t <-0, {f3,<j) <- (Po,a ) 
save (t, (3, a) 
while t < T max do 
w b ~ Af(0,l), w s ~ Af(0, 1) 

/3' <- /3 + ( M (<t) - D) (3 At + .^t^A + »JLJKtw- 

a'^a+ (-kfi(a) (3 + D (a* - a)) At + y^4§^ + ££■ + ff VStto" 
/3 -s— + 7. is an absorbing state for /3 
cr -s — |cr' — er min | + cr min % reflection at er min = — ^ s ln for a 
t^t + At 
save (t, f3, a) 
end while 



Algorithm 3: Normal approximation. 



The normal approximation (16) can be slightly modified in order to take into account the 
qualitative behavior of the SDE (17) near the axes. We propose the following scheme: 



A„ +1 = ft„ + [Mcr t J - (1 A K A & J D] f3 tn At 



Ki 



(lAKj p tn ) 2 D f) tn b 
K 4 w n 



*t» + [ ~ (1 A X 2 d tn ) k n(a tn ) (3 tn + D s- - (1 A if 5 <J tn ) D s] At 
At 



(lAK 2 a tn )2 fc^(s) ft,, 
K 2 



#3 



A' 5 



(25a) 



(25b) 



Indeed as (3 = is an absorbing state for the component f3 t of the SDE, instead of the standard 
Euler-Maruyama (16a), we can use (25a) where [•]+ is the positive part and w„ are i.i.d. Af(0, 1) 
random variables. 

Also, to take into account that the component a is reflected in er min = — |^ s' n we use the scheme 
(25b) where w a n are i.i.d. A/"(Q, 1) random variables. This discretization scheme was proposed in 
[4] in the context of the CIR diffusion process. In order to get a positive substrate concentration 
we can consider cr^ = a tn V or let <j min = in (25b). The simulation procedure is presented in 
Algorithm 3. 

Remark 4.1 (Scales and hybrid simulation) The three algorithms proposed here are valid at 
different scales. In the Gillespie algorithm all the detailed microscopic jumps of the dynamics are 
simulated. 
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cases 










K 2 


K 3 


K 4 


K 5 


o 


Hpf pfTTI 1 Tl 1 1 P 
















1 


"standard" 


CCLSG 


1.1 


10 4 


10 B 


10 b 


10 4 


10 b 




(see Figure 8) 


case 


1.2 


10 5 


10 7 


10 7 


10'"' 


10 7 








1.3 


10 7 


10 9 


10 9 


10 7 


10 9 


2 


"unstirred inflow/outflows" 


case 


2.1 


10 b 


10 b 


10 4 


10 4 


10 4 




(see Figure 9) 


case 


2.2 


10 7 


10 7 


10 5 


10 5 


10 s 






case 


2.3 


10 9 


10 9 


10 7 


10 7 


10 7 


3 


"fluid substrate" 


case 


3.1 


10 b 


oo 


oo 


10 4 


oo 




(see Figure 10) 


case 


3.2 


10 7 


oo 


00 


10"' 


oo 






case 


3.3 


10 9 


00 


oo 


10 7 


oo 


4 


"biological only" 


case 


4.1 


10 b 


10 4 


oo 


00 


oo 




(see Figure 11) 


case 


12 


10 7 


10 s 


oo 


00 


oo 






case 


4.3 


10 9 


10 7 


oc 


00 


oo 



(here "oo = 10 20 ") 



Table 3: Simulation cases. 

The idea of the Poisson approximation is to consider a time step A that should be small enough 
so that the different event rates barely evolve in the time interval [t,t + At], but large enough for 
the approximation to be worthwhile. Starting in x at t, the time step At is given by (24) but if it 
is less than a few multiples ofl/\(x) then the Gillespie algorithm should be preferred. 

Now the Poisson variables Vi of parameter Xi (x) At could be approximated by normal variables 
N(\(x) At, Xi(x) At) as soon as X l (x) At > 20. 

The simulation method can automatically switch from one algorithm to another one according 
to the scale. We can also imagine that different components of the state vector are simulated with 
different algorithms. 

5 Simulation study 

We present simulation results of the discretized diffusion model (25) with Monod and Haldane 
specific growth rates (2). The ODE (1) is integrated with a Runge-Kutta 1 scheme but the Eulcr 
scheme, corresponding to (25) with Ki = oo, gives very close results. 

In addition to the deterministic case (case with Ki = oo for all i), we consider 3 basic cases 
(see Table 3): 

"Standard" scales: -^2.3,5 = 100 x K\ t i corresponds to the "standard" case where the substrate 
concentration dynamics is closer to the deterministic case than the biomass concentration 
dynamics. 

"Unstirred inflow/outflows" scales: K\p, = 100 x #3,4,5 corresponds to the case where inflow 
and outflows are unstirred. 

"Fluid substrate" scales: #"2,3,4 = 00, in this case the substrate equation (17b) is deterministic, 
i.e. the substrate dynamics is in fluid limit. 

"Biological only" scales: K.3,,4, 5 = 00, in this case we consider that the randomness is only due 
to biological aspects of the system. 



lr The routine ode45 of Matlab, an explicit Runge-Kutta (4,5) formula. 
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Monod model 





set 1 


set 2 




k 


10 


10 


stoichiometric constant 


/Umax 


3 


0.5 


maximal growth rate (ft -1 ) 


D 


0.12 


0.4 


dilution rate (fa -1 ) 




0.5 


10 


input concentration (g/l) 


k~ : . 


6 


1 


half saturation constant (g/l) 




set 1 set 2 



Table 4: Sets of simulation parameters and the corresponding Monod growth rate functions s — > 
/i(s). The horizontal doted line is the maximum capacity fJ, max and the vertical doted line the 
asymptotic substrate concentration of the ODE corresponding to the non-washout case. 

5.1 A first comparison of trajectories 

We consider the set of parameters of Table 4 (set 1) for the Monod case in the "standard" scales 

Kx,a = 10 6 and #2,3,5 = 1q8 - 

In Figure 3 we present a simulation of the pure jump process with the Gillespie method. As 
expected, most of the events corresponds to small jumps of the substrate concentration. Before 
addressing the question of reliability of these algorithms, see next subsection, we first focus on the 
qualitative nature of the trajectories proposed by the various methods. 

In Figure 4, a simulation in a short time horizon of 0.1 (h) is proposed with Gillespie method 
(exact simulation) and the Poisson approximation method (tau-leaping) with a very small error 
control parameter e = 10 -6 . The corresponding trajectories are very similar though 10 6 events are 
needed for the Gillespie method and only 3500 time steps are needed for the tau-leap method. 

Figure 5 present a simulation on a realistic length of time 100 (h). Here only the Poisson 
approximation and the normal approximation are reliable. For the Poisson approximation we use 
e = 10 -3 (so that At is between 0.016612 and 0.034232, for 3004 time steps) and for the normal 
approximation we use 2000 (corresponding to At = 0.025). Again, the associated trajectories are 
very similar. 

5.2 Law of the concentrations at a given time t 

In Figure 6, we propose a Monte Carlo simulation to approximate the marginal densities of the 
biomass concentration B t and of the substrate concentration St at a given time t. We consider the 
set of parameters of Table 4 (set 1) for the Monod case in the "standard" scales #1,4 = 10 5 and 
#2,3,5 = 10 7 . 

We compute (S ( t 3 \ B ( t 3) ) for t = 3 (h) for j = 1 • • • 20000 independent Monte Carlo trials of the 
pure jump process (with the Gillespie method), with the Poisson approximation (tau-leap method) 
and with the normal approximation. For the tau-leap method we choose a constant time step. For 
"Poisson 1" and "Normal 1" we use a step of 0.05, for "Poisson 2" and "Normal 2" we use a step 
of 0.5. For each test, we compute the approximate density of S t and B t from the sample with a 
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0.026 

„ 0.026 

— 

': 

^ 0.026 
o 

| 0.026 
c 

8 0.026 
c 

| 0.026 
E 0.026 
■° 0.026 
0.026 



deterministic 
Jump process " 












0.04 0.06 
time (h) 



0.2602 

§ 0.2601 

f 0.26 

§ 0.2599 
8 

% 0.2598 

B 

= 0.2597 
0.2596 



deterministic 
- Jump process 



0.04 0.06 
time (h) 



Figure 3: Simulation of the pure jump process with the Gillespie Algorithm 1 for ^ = 10 6 
and 1^2,3,5 — 10 s and 1237928 events (blue: a realization of the jump process (Xt)o<t<o.i> green: 
the ODE (x(t))o<t<o.i) — the substrate dynamics (RIGHT) presents many small size jumps, the 
biomass dynamics (LEFT) less jumps but with higher amplitude. The corresponding phase-portrait 
is plotted in Figure 4 (LEFT). 




deterministic 
Poisson approximation 




0.2596 0.25960.25970.25970.2598 0.25990.2599 0.26 0.26 
substrate concentration (g/l) 



0.25960.25960.25970.25970.25980.25990.2599 0.26 0.26 
substrate concentration (g/l) 



Figure 4: LEFT : Simulation of the pure jump process with the Gillespie Algorithm 1 for 2£i,4 = 10 6 
and K2.3,5 — 10 8 and 1237928 events (blue: a realization of the jump process (Xt)o<t<o.i> green: 
the ODE (x(t))o<t<o.i) — the substrate dynamics presents many small size jumps, the biomass 
dynamics less jumps but with higher amplitude — the final time of simulation is 0.1 (h). RIGHT: 
same simulation with the Poisson approximation (^"t„)o<t„<o.i (tau-leap) with 3456 events for 
e = 1(T 6 (the time step is ~ 1.89 x 1Q" 5 J. 
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Figure 5: LEFT : Simulation of the Poisson approximation with the r-leap Algorithm 1 for = 
10 6 and 2^2,3,5 = 10 8 and 3004 events (blue: a realization of the jump process (-Xt„)o<t n <iOOj 
green: the ODE (x(t))o<t<iooA £ = 10 -3 , At between 0.028 and 0.0344. RIGHT: same simulation 
with the normal approximation (£t„)o<t„<ioo (diffusion approximation) with 2000 time steps with 
At = 0.05. 



kernel method. Hence we compare 5 probability density functions for each component St and B t . 
We also compute the empirical mean and standard deviation associated with the sample obtained 
from the pure jump process and we plot the associated normal density. 

Initial conditions are -Bo = 0.026 and Sq = 0.26, corresponding to the case of Figure 5, which 
is quite far from the equilibrium state. 

The conclusions are: 

• The two approximations (Poisson and normal) are very close to the exact simulation of the 
pure jump process; the approximation with a larger step 0.5 is slightly different. 

• The computation times 2 are: 

— for the exact simulation of the pure jump process: 5 h 45 min 32.6 s; 

— for the Poisson approximation: 33.2 s (with the time step 0.05) and 4.6 s (with the time 
step 0.5); 

— for the normal approximation: 0.7 s (with the time step 0.05) and 0.1 s (with the time 
step 0.5). 

In the present situation, where the parameters Ki are rather high, and for non-small con- 
centration of the biomass and the substrate, the exact simulation of the pure jump process 
( Gillespie method) should be avoided. 

• The resulting empirical densities are very close to normal densities and the solution of the 
ODE coincide with the mean of these normal densities. 

A second test is proposed in the case of the Haldane growth function (Set 2 of Table 5): see Figure 
7. The conclusions are the same as for the Monod case. 

2 CPU time on a 2.13GHz Intel Core 2 Duo with a RAM of 2 GB. 
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Figure 6: Empirical densities for the substrate and the biomass concentrations at time t = 3 
obtained with the exact simulation of the pure jump process (St,B t ) with the Gillespie method 
(blue line), with the Poisson approximation (St,B t ) with constant time step (red solid line for a 
step 0.05, red dash line for a step 0.5), with the normal approximation (at,/3t) with constant time 
step (green solid line for a step 0.05, green dash line for a step 0.5). The corresponding CPU time 
are respectively: 5 h 45 min 32.6 s, 33.2 s, 4.6 s, 0.7 s, 0.1 s. The simulation parameters are the 
Set 1 of Table 4 (Monod growth function) with "standard" scales ^ = 10 5 and ^2,3.5 = 10 . 
We compute the substrate and the biomass concentrations at t = 3 (h) and for j — 1 ■ ■ -20000 
independent Monte Carlo trials. The empirical densities are obtained with a kernel approximation 
procedure. We also compute the empirical mean and standard deviation from the sample of the 
pure jump process and plot with a thick grey line the corresponding normal densities: the match 
is very good. The vertical doted line the value of the substrate and biomass concentration at time 
t = 3 given by the ODE; again it matches the mean of all the empirical densities ( except the ones 
corresponding to the time step 0.5). 
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3500 



3000 



2500 



2000 



1500 



1000 



500 



substrate concentration p.d.f. at time 0.2 



normal pdf 

pure jump 

Poisson 1 

- Poisson 2 

— Normal 1 
Normal 2 



450 



400 



350 



300 - 



250 



biomass concentration p.d.f. at time 0.2 



i ■ 1 — ^-j 1 1 o 

0.9405 0.941 0.9415 0.942 0.9425 0.943 




200 



150 



100 - 



0.58 



0.585 



0.59 



0.595 



Figure 7: Empirical densities for the substrate and the biomass concentrations St and Bt at time 
t = 0.2 obtained with the exact simulation of the pure jump process (St,Bt) with the Gillespie 
method (blue line), with the Poisson approximation (St,B t ) with constant time step (red solid line 
for a step 0.005, red dash line for a step 0.01), with the normal approximation (a t ,/3 t ) with constant 
time step (green solid line for a step 0.005, green dash line for a step 0.01). The corresponding CPU 
time are respectively: 3 h 54 min 16.9 s, 13.5 s, 6.7 s, 0.1 s, 0.1 s. The simulation parameters are 
the Set 2 of Table 5 (Haldane growth function) with "standard" scales = 10 5 and if 2,3,5 = 10 7 - 
We compute the substrate and the biomass concentrations at t = 0.2 (h) and for j = !■■ -20000 
independent Monte Carlo trials. The empirical densities are obtained with a kernel approximation 
procedure. We also compute the empirical mean and standard deviation from the sample of the 
pure jump process and plot with a thick grey line the corresponding normal densities: the match 
is very good. The vertical doted line the value of the substrate and biomass concentration at time 
t = 0.2 given by the ODE; again it matches with the mean of all the empirical densities. The initial 
condition (Sq,Bq) is the stable equilibrium solution of the ODE (1) that does not correspond to the 
washout: hence these empirical densities could be considered as good approximations of the limit 
distribution of the stochastic chemostat. 
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Haldane model 





set 1 


set 2 




k 


0.1 


0.1 


stoichiometric constant 


/Umax 


2 


8 


maximal growth rate (fa -1 ) 


D 


0.4 


0.4 


dilution rate (fa -1 ) 




1 


1 


input concentration (g/l) 


k~:. 


4 


17 


half saturation constant (g/l) 


k/ 


1 


1 


saturation parameter ) 




set 1 set 2 



Table 5: Sets of simulation parameters and the corresponding Haldane growth rate functions s — > 
/J.(s). The horizontal doted line is the maximum capacity /i mam and the vertical doted line is the 
asymptotic substrate concentration of the ODE corresponding to the non-washout case. 

5.3 About the scales parameters 

As we have seen, for large populations, the diffusion approximation £ tn = (Pt nl Pt n ) given by (25) 
is very close to the reference pure jump model Xt = (Bt,St). So we now propose simulations of 
the diffusion approximation in the case of a Monod specific growth rate according to the scales 
scenarios of Table 3 

• Case 1 ("standard"): see Figures 8 and 12. 

• Case 2 ("unstirred infiowoutfiows") : see Figures 9 and 13. 

• Case 3 ("fluid substrate"): see Figures 10 and 14. 

• Case 4 ("biological only"): see Figures 11 and 15. 

Figures 8 to 11 represent a simulation of a single trajectory in the 3 levels of scale: cases m, 1 to 
m, 3 (for m = 1 • • • 4). Figures 12 to 15 represent the result of 10000 Monte Carlo trials in the 3 
levels of scale: cases m, 1 and m, 2 (for m = 1 • • • 4). We represent the mean trajectory and the 
empirical law of £y at final time T. 

We can conclude that, at this level of population and scale: 

• The stochasticity is negligible only in the Case 4 ("biological only") and at the highest scale 
level (cases m, 3). 

• The ODE solution x(t) matches the (empirical) mean of the stochastic process at these scales 
(as the stochastic process is solution of a nonlinear equation, there is no reason for the mean 
of the stochastic process to coincide with the solution of the deterministic equations). 

Equivalent results have been obtained for the Haldane case. 
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Figure 8: Diffusion approximation, Case 1, Table 3 ("standard case") / Simulation of (25) with 
Monod specific growth rate (2) with the parameters of Table 4 — Time evolution of the biomass 
concentration (top left), time evolution of the substrate concentration (top right), phase portrait 
biomass / substrate concentrations (bottom) according to 4 cases: case 0, case 1.1, case 1.2, case 
1.3 (see Table 3). Cases (deterministic) and 1.3 are identical. 
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Figure 9: Diffusion approximation, Case 2, Table 3 ("unstirred inflow and outflows") / Simulation 
of (25) with Monod specific growth rate (2) with the parameters of Table 4 — Time evolution of the 
biomass concentration (top left), time evolution of the substrate concentration (top right), phase 
portrait biomass / substrate concentrations (bottom) according to 4 cases: case 0, case 2.1, case 2.2, 
case 2.3 (see Table 3). Cases (deterministic) and 2.3 are identical. 
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Figure 10: Case 3, Table 3 ("substrate fluid limit case") / Simulation of (25) with Monod specific 
growth rate (2) with the parameters of Table 4 — Time evolution of the biomass concentration (top 
left), time evolution of the substrate concentration (top right), phase portrait biomass /substrate 
concentrations (bottom) according to 4 cases: case 0, case 3.1, case 3.2, case 3.3 (see Table 3). 
Cases (deterministic) and 3.3 are identical. 
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Figure 11: Case 4, Table 3 ("biological case") / Simulation of (25) with Monod specific growth 
rate (2) with the parameters of Table 4 — Time evolution of the biomass concentration (top left), 
time evolution of the substrate concentration (top right), phase portrait biomass /substrate concen- 
trations (bottom) according to 4 cases: case 0, case 4-1, case 4-2, case 4-3 (see Table 3). Cases 
(deterministic) and 4-3 are identical. 
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Figure 12: Cases 1.1 and 1.2 Table 3 / Sampling 10000 Monte Carlo trials of the law of (/3 tn , a tn ) 
for t n — 100 — The deterministic solution and the mean of the sampled trajectories coincide - 
The final law is represented by the sample and by the contour plot of the corresponding kernel 
approximation of the p.d.f. 




Figure 13: Cases 2.1 and 2.2 Table 3 / Sampling 10000 Monte Carlo trials of the law of (/3 t;i , <r tri ) 
for t n = 100 — The deterministic solution and the mean of the sampled trajectories coincide - 
The final law is represented by the sample and by the contour plot of the corresponding kernel 
approximation of the p.d.f. 
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Figure 14: Cases 3.1 and 3.2 Table 3 / Sampling 10000 Monte Carlo trials of the law of [fit n , &t n ) 
for t n = 100 — The deterministic solution and the mean of the sampled trajectories coincide - 
The final law is represented by the sample and by the contour plot of the corresponding kernel 
approximation of the p.d.f. 
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Figure 15: Cases 4-1 md Table 3 / Sampling 10000 Monte Carlo trials of the law of ' (0t„tO't„ 
for t n — 100 — The deterministic solution and the mean of the sampled trajectories coincide. 
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5.4 Comparison with the Imhof-Walcher model [15] 

We compare the processes £ t = (/3 t ,er t )* given by the diffusion approximation model (17) with 
the one given by the ad hoc model (21). The parameter are: K = 1, /i max = 1, D = 0.5, 
S"» = 8, fc 3 = 0.5, final time T = 20, At = 0.02, 20000 Monte Carlo trials; K, = 10 5 for (17) and 
c b = c s = 0.02 for (21). The parameters are chosen so that the biomass concentration evolves from 
0.5 to about 7.5, and the substrate concentration from 5 to about 0.5. Also the limit distribution is 
lesser than 1 in the substrate and greater than 1 in the biomass. Indeed one of the main difference 
between (17) and (21) is than for state values less than 1 (resp. more greater than 1) the noise 
variance for the first model is greater (resp. lesser) than the noise variance for the second model. 
This example illustrates clearly that the two models differ substantially. 

6 Discussion 

We started from a reference pure jump model X t , described by rates /jumps structure of Table 1 or 
as a solution of the stochastic differential equation (10). The martingale decomposition (41) clearly 
describes that the dynamics of X t is the combination of the classical deterministic dynamics of the 
chemostat (1) plus martingale terms with coefficients \j\fKl and with explicitly known quadratic 
variations, see (42). These quadratic variation terms allow us to assess the difference between the 
stochastic model and the deterministic one. 

We presented the explicit Monte Carlo simulation procedure, called Gillespie method, for the 
process X t . In standard cases, that is for high population levels (i.e. Ki large), this procedure 
is not feasible as it requires us to simulate too many events. In this case, we presented the 
Poisson approximation (25b) and the normal approximation (25), both in discrete-time. These 
approximations are valid only for large populations, i.e. about the axes, it is necessary to return to 
the pure jump process representation. In the application discussed here, the Poisson approximation 
is of little interest: it is more time-consuming than the diffusion approximation and valid only on 
a very limited scale range between the pure jump model and the normal approximation model. 

In contrast with previous stochastic chemostat models [26, 7, 15] where the stochasticity was 
introduced according to an ad hoc approach, in the present work we propose a family of models 
where the structure of the noise emerges from the very dynamics and where the scale parameters 
can be tuned according to the problem under interest. In particular it allows us to propose 
hybrid models where the cell population dynamics features stochasticity as the substrate is in fluid 
dynamics (ODE), corresponding to the Case 3 of Table 3. This kind of model has already been 
proposed in [11] in a three trophic levels case where the stochasticity appears only in the top level 
trophic as a stochastic logistic model and with fluid limit dynamics for the two other levels; it 
also has been proposed in [3] with a pure jump process for the biomass dynamics and a fluid limit 
for the substrate. This approach can also be related to coupled slow/fast reactions in stochastic 
chemical kinetics [13, 2]. 

The approach proposed here can be applied to any model of population dynamics especially in 
cases of difference of scale between the different dynamics (e.g. cell/substrate). The dynamics of 
interacting populations cannot be modeled by a single model but rather by a family of models whose 
domain of validity depends on the scale at which the dynamics are considered. For example the 
normal approximation model represented as stochastic differential equations (17) or the ODE model 
(1) are valid in high population levels, hence using such models to infer extinction characteristics 
like extinction time and extinction probabilities is not valid. This was already noticed by [22] and 
[27]. 

In most standard population scales of the chemostat the ODE model is justified. Also, the ODE 
framework proposes analysis, control and optimization tools that are more accessible than the one 
of the SDE context. Though, as seen, the stochasticity cannot be neglected in many situations. 
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Figure 16: Comparison of the processes £t = (/3t,<Jt)* given by the diffusion approximation model 
(17) and by the ad hoc model (21). Evolution of the biomass concentration (top) and substrate 
concentration (center) during time; final joint density approximation of concentration (bottom). 
The two models differ substantially: compared to the diffusion approximation, the ad hoc model 
overestimates (resp underestimates) the noise variance for concentration greater (resp. lesser) than 
1. 
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This stochasticity could be of small intensity in the present single species/single substrate situation 
but could deeply perturb multiple species/multiple substrates situations. The SDE model could be 
simulated at a small extra computational cost and offers a more realistic prediction tool. Indeed, as 
it can account for the variability of the experiments, the simulation of the SDE offers the possibility 
to explore in depth the potentialities of the dynamical systems. 

The SDE model is also more adapted for the confrontation to the data as it allows us to build 
a statistical model and the associated likelihood function. One of the next important steps, that 
we will investigate in coming work, will be to propose an adapted statistical procedure to estimate 
the scale parameters Ki, and in a second step to estimate the parameters (D, s ln ...). In the future 
we will also investigate the long-term behavior of these models as well as their optimal command. 
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Appendices 

A Representation of the process X t 

Infinitesimal generator 

We consider the Markov process X t with infinitesimal generator: 

5 

Acj>(x) = [(/>(x + Vl (x)) - <j>{x)\ (26) 

i=l 

for all <f> : i-> K continuous with compact support [6, Th. 8-3.1]. The infinitesimal generator 
can also be understood in the following way: 

!Xi(x) At + o(At) , if y = v^x) for i = 1 • • • 5, 

1 ~ Ei=i \ (x) At + o(At) , if y = 0, 
o(Ai) , otherwise. 

or as 

t-+o t 

where Xf is the process X t starting from x. It can be rewritten as: 

A(f>{x) = X(x) I [<j>(y) - 4>{x)] p(x, dy) 

with 

A(x) d = f ^AiW, (27) 



i=l 
o 



p(x,dy)=y2Xi(x)8 x+1/i ( x) (dy) with Ai(a;) = . , (28) 

We define the jump times: To = and 

r„ = inf {t > r„_i ; X t 7^ } 

and the embedded jump chain 

v — Y 

It is well know that (i) Y n is a Markov chain on with transition probability p(x, dy); (ii) for all 
ti ^ 1, conditionally on Yo, • * • ? ^n— l; the holding times tx — To, ... , T n — T n _i are independent and 
exponentially distributed of intensity parameters X(Yq), . . . , A(i / rl _i), see [21]. These properties 
are at are the basis of the Gillespie simulation algorithm (see Algorithm 1). 

Non-explosion and existence of moments 

To study the non-explosion and existence of moments, we define the mean jump size function: 

r 5 - 

m K (x)= (y - x) p(x,dy) = V v l {x) Xi{x) . (29) 
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and, for p > 1 

r 5 

\m K \ P (x)= \v-x\*p(x,dy)=y2\u i (x)\*X i (x)- (30) 

Note that from (14): 

fxix) = A(x) mif (ar) (31) 
and if we replace Vi(x) by i/j in (31) we get: 



where f(x) is the right-hand-side function of the ODE (1). 

fx(x) is the instantaneous mean of the process Xt, more precisely we will show in Proposition 
A. 3, that K(Xt+At\Xt = x) ~ Jk{x) At. As proved in the next lemma /k{x) is essentially f(x), 
so locally in time the mean of the process X t behaves like x(t). 

Lemma A.l Consider fi(x) defined in Table 1, vi defined by (9) and fx (x) defined by (14). First 
\ u i( x )\ < Then let 

K K = {x = ( b s )eRl;b<^ors<^or S <^}. 
For x TZk, v%(x) = v i; m^(x) — m(x), fx (x) = f(x). For all x € M. 2 :: 

\vi(x)-Vi\<— — , i=l,...,5, (32) 

nmxj=i...5 Kj 

\m K {x) - m(x)\ < — — , (33) 

mmj = i...5 Kj 

\f K (x)-f(x)\ < (1- K 2 s)+kp(s) b+ {1-K 4 b)+Db+ {1-K 5 s)+Ds (34) 

so that Jk(x) — > f(x) when Ki — > oo for all i = 1, . . . , 5 and this convergence is uniform on any 
compact set of (0, oo) 2 . 

Proof For x ^ TZk it is clear that Vi(x) = vi so /k(^) = /(#)• For all x, we have v\{x) = v\ and 
v%{x) = V3 and: 





-V2 = 


K 2 1 






-K 2 s) 


v 4 (x) 




^1 


v (1 " 


-K 4 b) 



u 5 (x) 


-v$ = 








-Kb s) 



so we get (32). The following assertions of the lemma are straightforward. □ 

Non-explosion and existence of moments are given by the following result: 

Theorem A. 2 (Hamza and Klebaner [12]) Suppose that X(x) > for all x <G M.+ and that 
there exist C > and p > 1 such that E(|Ao| p ) < 00 and 

X(x)\f K \ p (x) < C(l + \x\ p ), VxeR 2 + (35) 

then the Markov process (X t )t>o is non- explosive, that is r„ — ¥ 00 a.s., and for all T > there 
exists C > s.t. E(\X t \P) < C~for all t < T. 
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Indeed X(x) = J2i M 1 ) > K 3 D > and 

5 
i=l 

H(s)b kp,(s)b Ds" 1 Db Ds 



< 



Kl~ x K?- 1 K\- x Af 1 

so that (35) is fulfilled for all p > 1 . Hence the Markov process X t with infinitesimal generator A 
defined by (26) is non-explosive and X t admits moments of all order for all t > 0. 

Representation for the process X t 

We now give a representation for the process X t as a solution of a stochastic differential equation 
driven by random Poisson measures: 

Proposition A. 3 The process X t is defined for all t > and it is solution of the jump SDE (10) 
where N l are independent random Poisson measures with intensity measure du x dv (Lebesgue 
measure). 

Proof We first verify that the stochastic integral in (10) is defined. According to [14, §11-3], this 
integral is defined if: 

/>£ r-OO 

E \vi(X u )\ l{v<x % (x u )} dvdu < oo . 

Jo Jo 

We have 

ft poo 



JO 



E / \L>i(X u )\ l {t) < A . (Xn)} dudu 



<-^e/ A 2 (A u )dw < C [ (l+E|A„|)du 
K i Jo Jo 

which is finite according to Theorem A. 2. □ 



Consider the centered random Poisson measure: 



N l (du x du) = N*(du x dv) — du x dv . 



According to [14, § II-3] 



f-t r-oo 

M t id = f / / Vi(X u -)l {v < Xi(x )} JV*(d«xd«), i = 1 5 (36) 

Jo Jo 

are five independent square-integrable martingales with finite moments of all orders. Let: 

M t ^Y.Ml. (37) 

i=l 

From (10): 

X t = X + f f K (X u ) du + M t . (38) 
Jo 
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We want to study the behavior of the martingales Ml as the Ki — » oo. First note that 



Ml 



( m o0 



for i = 1,4 and Ml 



( m 



for i = 2, 3, 5 with 



1 def 

to, = 



2 dot 

to, = — 



3 d « f 

TO? 



j />* />oo 

IT 1 {«<Ai(x„ _)} J /V 1 (du x du), 

A i Jo Jo 

If'/" 00 

— / (lA/f 2 V)l{»<A 3 (x _ )} iV 2 (d U x d«), 

A 2 Jo JO 

1 /"* Z" 00 
IT / / 1 {^<A3(x„_)}^ 3 (dwxdw), 

A 3 Jo JO 



4 dcf 



— / / {lhK A B u -)l {v < Xi{Xu _ )} N\Aux&v), 



K i Jo Jo 



5 d <= f 
TO = - 



1 



K, 



5 Jo Jo 

t poo 



(I A K 5 S u -)l {v < MXu _ )} N 5 (du x dv) . 



(39a) 
(39b) 
(39c) 
(39d) 
(39e) 



As m\ is of the form m\ — J J °° "f l (u , v) N l (du x di>) then the associated predictable quadratic 
variation is (m 4 ) t = J* / °° [y l (u, v )] 2 dvdu so we can easily check that 



1 /"* 

m 1 )* = -rr \ n{S u )B u du, 
Ki Jo 



m 2 )t 



AX 2 S' tl ) 2 fc/x(5*„) J B„dw, 



Ds ln t, 



in 



AK 4 B U ) 2 DB U du, 



m 5 ) t 



K 5 Jo 



(1AK 5 S U ) 2 DS u du. 



Define 



'/Q mi = \JK, 



i / 7*(u ,v) N^du x du) . 

J[0,t]x[0,oo) 



(40) 



So we obtained the following representation of the process X t that emphases the dependence on 
the Kr. 



dB t = (fx(S t ) B t - (1 A Ki B t ) D B t ) dt + -±= dm\ + -±= dm 
dSt = ( — (1 A K<i St) k fi(St) B t + D s" 1 - (1 A K 5 S t ) D S t ) dt 



(41a) 
(41b) 
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where m\ are independent square integrable martingales with the following quadratic variations: 



m 1 )^ 


/ (j,(S u )B u du, 
Jo 


(42a) 


m 2 ) t = 


f (lAK 2 S u ) 2 k[i(S u )B u du, 
Jo 


(42b) 


m 3 ) t = 


Ds ,n t, 


(42c) 


m )t = 


f (1AK 4 B U ) 2 DB u du, 
Jo 


(42d) 


m 5 ) t = 


[ (lAK 5 S u ) 2 DS u du. 
Jo 


(42e) 



B Existence and uniqueness for a solution of the SDE (17) 

To prove that the system (17) admits a strong solution and pathwise uniqueness holds we use the 
results of [24, p. 134] or [5]. Let V N = [jj,N] x s in + jj,N] and rewrite (17) as: 

dSt = f(Zt)dt + g(Z t )dW t (43) 

where £ t = ((3 t ,o- t )* and Wt = {W^Wf)*. The coefficient g is not Lipschitz continuous, it is 
globally Lipschitz on T>n for all N, but we can find Lipschitz continuous coefficients /at and gpj 
such that: 

In (0 = /(£) for £eV N , f N (O = for £ V 2N , 

»w(0 = 5(0 for £ G Pjv , ffjv(0 = for £ £ P 2A r . 

The the system 

d^ = / w (ef)dt+ 5J v(^)dw t 

admits a unique strong solution £ N . Let: 

T N = M{t > ; £f £ Pjv} . 
Hence if ./V > M then £ w = £ M for all t < Tm, so we can define a process £°° such that: 

and will be solution of (43) up to the explosion time Too = limjv-^oo T/y. Then by stability property 
and by the fact that the solution cannot cross the boundary of V we get = oo a.s. which proves 
the strong existence and pathwise uniqueness defined for all t > 0. 
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